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1 • INTRODUCTION . 

This paper addresses the rather narrow subject of find- 
ing airfoil and wing shapes that are free from shock waves 
even though the local flow speed exceeds the speed of sound. 

This research stems directly from an invitation we made to 
Helmut Sobisiczky of the Deutsche Forschungs- und Versuchsan- 
stalt fur Luft- und Raumfahrt in Gottingen, West Germany, to 
spend the academic year 1977-1978 at the University of 
Arizona. My specific proposal to Sobieczky was that we 
collaborate on replacing his complicated analog computations 
of solutions to the hodograph equations by a fast elliptic 
solver in order to generate shock-free airfoil designs more 
effectively. The first part of this paper addresses this 
study. The second part addresses a much more efficient pro- 
cedure that we now use to the same end. This second method 
is a result of Sobieczky 's brilliant idea of a fictitious gas 
for finding shock-free airfoils directly in the physical 
plane. 

The aerodynamic efficiency of turbojet-powered aircraft 
is proportional to the flight Mach number times the lift-to- 
drag ratio, viz., M^.L/D. In addition, the amount of return 
an aircraft provides on the investment in it is proportional 
to the flight speed, and hence, to the flight Mach number. 

At supersonic Mach numbers shock waves are always present and 

give rise to a wave drag that adversely affects the aero- , 

A 



dynamic efficiency, and they are also the cause of the sonic 
boom. As Albert George and I showed some time ago [17], the 
sonic boom has an irreducible minimum for a given aircraft 
weight and length, and this is large enough to preclude com- 
mercial operation of supersonic aircraft over populated areas. 

The rapid rise of petroleum prices has driven the cost 
of fuel rapidly upward and it now comprises more than 50 per- 
cent of the direct operating cost of an aircraft. If air 
transportation is to remain affordable, we must draw on a 
number of technologies to provide a substantial improvement 
in the efficiency of transport aircraft. Fortunately, there 
are a number of technologies that, when combined, should 
provide a 40 percent improvement in seat-miles per gallon 
[18]. For a single DC-10 size aircraft, each 1 percent 
improvement in M^L/D saves about $100,000 per year in fuel 
costs at present prices [9]. Indeed, by the year 2000, we 
should achieve 90 seat-miles per gallon for intermediate 
range flights. 

One, albeit small, element of this improvement is better 
aerodynamic efficiency through special airfoil and wing 
designs that allow flight at supercritical Mach numbers (i.e., 
subsonic flight Mach number high enough that local regions of 
the flow are supersonic) without shock-waves. This avoids 
the drag associated with this shock wave and, more impor- 
tantly, the boundary-layer separation that occurs if such 
shock waves become very strong. 

This goal was once thought unattainable. In the mid- 
1950 's Morawetz [11] showed that smooth, i.e., continuous, 
flows with embedded supersonic regions were mathematically 
isolated from one another. That is, any small change in the 
Mach number, or any arbitrary change in the airfoil surface 
embedded in the supersonic region, would result in the 
appearance of a shock wave in the flow. This led a number of 
people, including the author, to conclude that such flows 
would be of no practical interest. Others were more wise. 

Wind tunnel studies by Pearcey [14], and by Whitcomb [26] 
showed that such flows did exist and were of practical im- 
portance. Subsequent theoretical studies by Nieuwland [13] and 
Bauer, Garabedian and Korn [1] led to techniques for finding 



shock-free airfoil shapes. Boerstoel [2] later improved this 
technique. Sobieczky [ 19 ] developed yet another method that 
employed a rheo-electric analog computer to solve hodograph- 
like equations . These earlier contributions are not reviewed 
here . 

The first part of this paper deals with the inverse 
problem: given an airfoil pressure distribution, find the 
airfoil that has that pressure distribution. And it also 
deals with the extension of this procedure to airfoils with 
an embedded supersonic region. This we call indirect design. 
The results are a direct consequence of Sobieczky 's rheo- 
electric analog design method where the emphasis was on find- 
ing shock-free airfoils. This extension was conceived and 
used by Sobieczky [20] and further developed by Hassan [5]. 

The second part of this paper treats Sobieczky 's gener- 
alization of this process [21]; it provides a tool for the 
routine design of shock-free airfoils and wings. 

I I . INDIRECT DESIGN OF AIRFOILS . 

The basic problem we address here is to find an airfoil 
in two-dimensional, irrotational , flow that has a prescribed 
pressure distribution. This flow is described by a stream 
function, if), that guarantees the conservation of mass,, that 
is, 

curl(ipk) = pq 
so that 

div(pq) = 0, 

where 

• P/P*, = [1 + ^M^l - q 2 )] 1/Y_1 , 

and the symbols have their usual meanings. Because the flow 
is irrotational there is a velocity potential, <f> , such that 

q = grad <p . 

In the direct problem we seek solutions of 

= 0 ( 1 ) 

where 


n 


0 



on the airfoil and <j> is asymptotic to the potential for a 
free stream flow past a vortex, of strength equal to the jump 
in <p at the airfoil's trailing edge, located inside the air- 
foil. Or we may solve a similar equation, viz., 

Q 2 i> = 0 (2) 

with 


i|> = 0 

on the airfoil and with similar asymptotic boundary condi- 
tions. Here and are simple quasi-linear operators of 
second order. If we interchange the dependent and indepen- 
dent variables in (1) and (2) , then we will have linear 
equations for the coordinates. It is much simpler to intro- 
duce the hodograph variables q, the flow speed, and 0, the 
flow deflection angle as independent variables. This leads 
to (see, e.g. , [25]) 


L^(f> = 0 or = 0, (3) 

where and L 2 are linear, second-order operators. But we 
do not know the airfoil boundary in the q ,0-plane. However, 
in our approach we wish to specify the pressure, or equiva- 
lently q, and find the airfoil. Figure la is a sketch of a 
supercritical flow past an airfoil that has lift. In this 
case the hodograph image of the flow will be two-sheeted and 
appear like that sketched in Figure lb. Actually, it is 
mathematically convenient to introduce the variable 


v = 



dq 

q 


in place of q in equations (3) to reduce equations (3) to 
their canonical form. For the stream function we have 


^vv ± ^00 K ( v ) °' 

where ± is appropriate if M £ 1, and 

K(v) = J\ 1 - M 2 |/p. 


The subscript ( ) * refers to the sonic state where q = a = 
a* and M = 1. The variable v is the usual Prandtl-Meyer 
function for M > 1. Sobieczky refers to these variables, 
viz., v, 0, as the rheograph variables because equation (4) 



a\oj 



can be solved by electrostatic means if one supplies a con- 
ductor of varying conductivity, as he did by punching holes 
of varying size in conducting paper. But the general form of 
equation (4) is invariant under a conformal mapping and we 
can clearly map the subsonic region of the flow (see Figure 
lb) in the q,9 or v,0 planes into a unit circle. This simple 
device, suggested by Sobieczky [ 2 0 ] , and followed by Hassan 
[5], provides us with a simple boundary value problem that 
may be solved efficiently with a fast direct solver. Thus, 
we map the subsonic portion of Figure lb into the unit circle 
of Figure 2 with an unknown mapping of the form 

S = re iaJ = f (v + i0 ) . 

We insist that this mapping be conformal and hence v(r,to) and 
0(r,oj) are harmonic: 

2 2 

V v = V 0 = 0. 

If we give the pressure on the mapped airfoil where the flow 
is subsonic, then we know v there. And on the sonic line 
v = 0 .. So v is known on the unit circle and we can find v 
and 9 to within an arbitrary constant by Fourier series . We 
must be careful to take into account the behavior of v as 



q -* 0, viz., v ■= log q. Given the flow conditions at infin- 
ity, we know v^, 0 , and this locates the point at infinity, 
I. The stream function is singular there, and hence we use 
the bilinear transformation to map I to the origin of the a 
plane 


Now on the airfoil boundary ip = 0 and we may provide a guess 
for if) (w ' ) on the sonic line. Thus, we iteratively solve 


V 2 t|; n 


= [■ 
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with as e -*■ 0 and tf/(l,co') known, using a fast Poisson 

solver. Hassan [5] employs the sixth-order solver of Roache 

-1/3 

[16]. Because on the sonic line ifi , ( 1 , to ' ) « g (w ' ) ( 1 - r') , 

Hassan uses the isotach for = 0.99 in place of the sonic 

line and finds the sonic values by Taylor series. The 
physical plane image of the airfoil and sonic line combina- 
tion is simply found by integration around the unit circle 
using (see, e.g., [16]) 
ie(co-) 

dz = [<J> w , ( 1 ,( 0 ' ) - iip r , (l, ( o>)]dco'. 


This does not necessarily provide a closed airfoil and the 
input must be adjusted to obtain a satisfactory solution [6]. 
Now the given pressure distribution is associated with known 
points on the airfoil and the input pressure distribution may 



Figure 2. The (r,u) plane image of the subsonic portion of 
the airfoil surface and the sonic line. See 
Figure 1. 



have to be adjusted to achieve a satisfactory pressure dis- 
tribution on the airfoil. 

For supercritical flows the data on the sonic line is 
integrated using the method of characteristics to find the 
airfoil shape consistent with the sonic line data. To be 
specific , 

d<}> + K(v)dtj; = 0 (5a) 

on 

v + 0 = const, (5b) 

and the sonic line data $(0) and ip(0) are integrated to find 
the locus ip = 0, that is, the airfoil. Now the Jacobian of 
the map back to the physical plane may vanish before the line 
ty(v,9) = 0 is reached because M > 1, in which case a limit 




Figure 3. Comparison of the design pressure and sonic line 
with that computed for the designed airfoil. 



line [25] may appear in the physical plane. This means that 
the sonic line data is not compatible with shock-free flow, 
and new input data must be given. 

One typical result of the process is shown in Figure 3 
(from [20]). While this method is effective in finding 
shock-free airfoils, a more direct procedure for doing so is 
described in the next section. The main virtues of this pro- 
cedure are that it is an effective tool for finding sub- 
critical airfoils with desirable pressure distributions, and 
it led Sobieczky to his clever idea [21, 22] for the nearly 
direct design of shock-free airfoils and, more importantly, 
wings „ 


111 . SHOCK-FREE AIRFOIL DESIGN . 

Recall that in the previous section we used a fast 
Poisson solver iteratively to solve for the subsonic portion 
of the flow past an airfoil. For supercritical Mach numbers' 
this provided data on the sonic line and the supersonic flow 
was computed using the method of characteristics. If the 
sonic line data was inconsistent with shock-free supersonic 
flow, then the input was modified. Sobieczky [21] observed 
that this process could be mimicked in the physical plane by 
changing the equation relating the density to the flow speed 
when the flow was supersonic in such a way that the governing 
equation remained elliptic. Thus we solve 

div(pV<j>) = 0 (6a) 


where 


P/P 

' ^ 00 


[1 + 


Y-l 

2 


M*(l - 




q <_ a 


(6b) 


P f /Pco, q > a. 

If we wish equation (5) to be elliptic, and we limit our 
attention to fictitious densities of the form p^ = p^ (q) , 
then we must have 


d.(p f q) 

dq~ 


> 


0 . 


Thus we compute the flow past a baseline airfoil which has 
been selected on the basis of its subcritical performance. We 
may wish to alter it slightly in advance of our calculation, 



using equation (6) to give it more thickness to compensate 
for a subsequent reduction in thickness. This fictitious 
solution satisfies the boundary conditions and the correct 
equations where q < a. We take the data, viz., $ and 6, that 
the solution provides on the sonic line and generate ip ( 0 ) 
from the known relations between <f> and if). We then recalculate 

the supersonic domain using equation (5) . Of course, if a 
limit line intervenes we must alter our baseline airfoil, the 
free stream Mach number, the angle of attack, or the ficti- 
tious density relation (6b) and repeat the process. 

Figure 4a compares the pressure on an NACA 64A410 air- 
foil before and after it has been subjected to this design 
process; it also shows the sonic line (and shock) on the 
original 64A410 computed using the correct density law, and 
the sonic line computed using the fictitious density law (6b). 
The calculation of the supersonic flow by characteristics 
with the correct density-flow speed relation defines a new 
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CL 

0.7029 
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CO 
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Comparison of Ihe pressure coefficients and sonic lines for Ihe 
baseline NACA 64A410 and Ihe shock-free airfoil obtained from it by 
Ihe direct design procedure. 



HACH = .720 ALPHA : 0.40 
OESIGN ( 0 I ANALYSIS ( * I 
CL 0.7012 0-7029 
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CM --1395 -.1397 

Comparison of Ihe pressure coefficient and Ihe sonic line 
obtained by Ihe design calculation that modifies the airfoil shape with 
those obtained by computing Ihe flow past Ihe modified airfoil. 


Figure 4. Airfoil designed using the fictitious gas p f = p*. 



airfoil shape in the supersonic zone and predicts the in- 
dicated pressure on this new shape. A recalculation of the 
flow past the resulting airfoil is compared with the pre- 
dicted pressures in Figure 4b. The algorithm used is 
Jameson's FL06 [7] which has proved to be very reliable for 
inviscid flows. Of course the lift coefficient and thickness 
of the designed airfoil are less than those of the original 
airfoil. We may find additional designs by fixing the lift 
coefficient and varying the Mach number, by fixing the Mach 
number and varying the angle of attack, or by changing the 
gas law. Frequently we have used 
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V SCALE = S«tX SCALE) 

Shock-free airfoil shapes for fixed lift coefficient C, = 0.70 
and varying Mach number. The fictitious gas has a constant density in 
the supersonic domain (/> = 0). The baseline airfoil is an NACA 
64A4I0. 



1 SCALE = 5*1 X SCALE! 

Shock-free airfoil shapes for fixed Mach number =0.72 
and varying lift coefficient. The fictitious gas has a constant density in 
the supersonic domain (f > = 0). The baseline airfoil is an NACA 
64A410. 

Figure 5. Shock-free airfoil 


(7) 



v SCALE = 5«l X SCALE I 

Shock-free airfoils for fixed Mach number M „ = 0.72 and lift 
coefficient C , =0.70, varying the exponent P of Eq. (Ic) and thus 
changing the density's dependence on flow speed. The baseline airfoil 
is an NACA 64A410. 
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Parameter space explored for the shock-free airfoils that can 
be obtained when the baseline configuration is an NACA 64A410 
airfoil. 


shapes for varying , C , and P. 



where P < 1. The results of such a study are shown in Figure 
5. They indicate the wealth of solutions available by this 
technique. It also depicts shock-free designs found in P, 

C ,M space . 

Li 00 

Ilia. VISCOUS EFFECTS . 

Real flows are of course viscous and even for very 
large Reynolds numbers we must take the boundary layer dis- 
placement effect into account. In the main, such effects are 
weak and we may calculate the inviscid flow and use its pres- 
sure to compute the boundary displacement thickness, 6*. 



Figure 6 . 


Comparison of C 's and sonic lines for two 11.7% 
thick airfoils.^ Airfoil (a) is a thinned GA(W)-2; 
airfoil (b) is a shock-free design. 



There are local, strong inviscid-viscous interactions at the 
foot of a shock, at the trailing edge, and when the boundary 
layer separates from the airfoil, and in these regions we 
must couple the computation of the inviscid flow field with 
that of the boundary layer. The strong interaction that 
occurs at the trailing edge has been coupled with Green's lag 
entrainment model of the integral boundary layer equations by 
Melnik and his co-workers at Grumman [10] in the Grumfoil 
algorithm. If we make the not so brash assumption that the 
boundary layer displacement thickness for a shock-free design 
will grow in much the same way as it does for the fictitious 
gas we may try to find shock-free designs in the presence of 
viscous-inviscid interactions, even when they are strong. 

Nebeck [12] has conducted such a study with the Grumfoil 
algorithm and markedly improved the design of the Va-2 air- 
foil designed by the German aerospace industry. A subsequent 
study by Cosentino (private communication) using the GA(W)-2 
airfoil provided an 11.7 percent thick airfoil that he 
estimates to have 73 less counts of drag than an 11.7 percent 
thick version of the GA(W)-2. His airfoil and the pressure 
distribution computed using Grumfoil are compared with the 
thinned GA(W)-2 and its pressure distribution in Figure 6. 

IV. SHOCK-FREE WING DESIGN . 

When this airfoil technology is applied to wings, 
especially swept wings, the results are not very satisfactory. 
Indeed, typical results are sketched in Figure 7, which shows 
that three-dimensional effects will lead to shock waves 
unless they are explicitly taken into account in the design 
process. Conceptually, at least, we can imagine the process 
of Section II extended to three dimensions where we tackle a 
problems like that sketched in Figure 8a. But we must recall 
the sketch of Figure 8b. In two dimensions we cannot dis- 
tinguish the forward facing Mach lines, PA, PC, emmanating 
from point P, from the upward facing Mach lines PC, PB that 
also originate at P. But once we consider the flow to be 
three-dimensional, we surely can distinguish the fore Mach 

\ 

cone PAA'C'C from the surface PCC ' A' ADD 1 B ' BP . In other words, 
while we cannot distinguish the time-like direction for 
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yy 
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Figure 7. Three-dimensional effects on a wing constructed 
from shock-free airfoils. 



surface that results in elliptic equations. Sketch of Mach conoids for a two-dimensional flow. 

Figure 8. Ill-posed boundary-value problem that arises for 
wings . 


clearly only x is time-like in 


-<j> + <f> + <p = 0. (8) 

When we extend Sobieczky's fictitious gas concept to three 
dimensions we are essentially trying to solve equation (8) by 
giving data on the x,y-plane as sketched in Figure 9, and 
trying to find the solution for negative z. This problem is ' 
clearly ill-posed. The modal solution to equation (8) is 


Mach cones of the linear wave equalion and their intersection 
with the plane where the initial values are given. 

Figure 9. Boundary-value problem for equation (8). 

exp[i(k.x + k^y + k^ - k^ z) ] 

and if the wave numbers of interest in the y direction, k^, 
exceed those of interest in the x direction, k^, the solution 
will grow exponentially. Indeed, L. E. Payne (private com- 
munication) has shown that for any initial data the solution 
will eventually grow exponentially in z. But for most prac- 
tical solutions the wave numbers of interest in the x, or 
chordwise, direction are much higher than those in the y, or 

spanwise, direction. One can, of course, precondition the 
initial data to suppress the instability associated with this 
ill-posed problem, but we have not found this to be necessary. 

One of the frequently used algorithms for three-dimen- 
sional wings is FL022, which was developed by Jameson and 
Caughey [8], It is a nonconservative code and will not cap- 
ture shock waves correctly. But for our purpose of shock- 
free wing design. it is perfectly adequate. Here we must 
introduce a fictitious sound speed rather than a fictitious 
density. Such sound speed laws and their relation to the 
fictitious density laws may be found in [4]. We compute the 
supersonic flow by considering the conservation of mass, 

(a - u ) u + ( a - v ) v + ( a - w w 
x y y 

- 2uvv - 2uwv = 0 
X z 


and two of the three xrrotationalxty condxtxons, vxz. , 

w - u =0 and w - v = 0 , 
x z y z 

as a system of three equations for the vector U = transpose 
(u, v, w) , viz . , 

AU + BU + CU =0. 

—x — y — z 
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Figure 10. Three-dimensional shock-free flow over a rec- 
tangular wing of aspect ratio 10 at = 0.70. 
The NACA 64A410 airfoil section was used as a 
baseline for the wing. 


We then spline fit the initial data to compute and and 

find U to calculate the data of the next z-level. In 
— : z ^ 

general, then, at the k -level 


U, = U, , - [C -1 AU + C'" X BU ], 

— k — k-1 L —x — y J k- 


- 1 , 


1/2 


AZ, 


-1 


where [ . . . ]r-i/ 2 i n ^icates a suitable average. Of course C 
may not exist, indicating the intervention of a limit surface 
in the supersonic domain. We continue until we find U on the 
wing surface. For small changes in wing shape a linear, 
first-order partial differential equation can be solved to 
find the new wing shape. Further details may be found in [3, 
4, 15 and 22]. Yu [27] has also applied this method to wing 
design using a finite volume code. 

Figure 10 shows the results of this design process for a 
simple rectangular wing with a half-span to chord ratio of 5. 
The basic wing section was an NACA 64A410. The predicted 
pressures on the original wing and those on the wing designed 
to be shock free are the result of the wing section changes 
indicated below the wing. 


IVa . VISCOUS EFFECTS . 

We have been able to include viscous effects through 
the Pablim algorithm of C. L. Street (private communication) , 
which couples the three-dimensional boundary layer code of 
Stock [24] with FL022. When we compute the flow past the 
original wing, we use the displacement thickness computed 
using the shock-free design. This is done to suppress the 
separation of the boundary layer and favors the original air- 
foil. We find that the new design, at least at its design 
point, has an "aerodynamic efficiency" or range factor that 
typically exceeds a non-shock-free design of the same thick- 
ness by 5-15 percent. Such gains are of great practical 
importance and achieved with minimum effort. A study we have 
made on improvements to the Learjet Century series aircraft 
is given in Ref. 3. The GA(W)-2 airfoil used here is 
probably already superior to the present Learjet airfoil 
section. Certainly the redesigned wing provides an improve- 
ment in supercritical performance and, because it is based on 
the GA(W)-2 airfoil, should have good low speed performance 
as well. 



V. THE DESIGN PROCESS . 

The design process for shock-free wings must begin with 
an airfoil and wing planform that provides good low speed and 
subcritical performance. The art of the designer is in 
choosing a baseline airfoil and fictitious gas law that will 
maximize M^L/D without compromising the established subsonic 
performance. Here experience is important, as the designer 
must know what changes in the baseline airfoil and fictitious 
gas law will produce the desired effects. But current compu- 
tational tools provide this knowledge at little computational 
expense. Indeed, an undergraduate student is conducting the 
investigation of an improved wing for the Lear jet Century 
series aircraft. 

VI. CONCLUSIONS . 

While the flows we seek represent mathematically iso- 
lated solutions, there is a great wealth of such solutions 
and they are easy, and inexpensive, to find. Spee [23] has 
shown that In two dimensions if the local Mach number 
exceeds 2/^(3 - y) =1.58 then the flow must be unstable to 
small unsteady disturbances. Perhaps such shock-free solu- 
tions exist with higher local Mach numbers, but we have not 
found them yet and the condition M =1.58 could be a goal the 
designer should try to achieve. We don't know yet how the 
drag depends on the flight Mach number at a shock-free design 
point, nor have we established any theoretical limits on the 
Mach number, lift coefficient and thickness, beyond which 
shock-free designs cannot be found. Such theoretical limits 
could be most useful to the designer. 

Let me note in closing that while the ideas fundamental 
to this investigation are simple, they are mathematically 
motivated, and they draw upon a wide body of traditional 
applied mathematics. Also, their implementation depends on 
computer algorithms carefully constructed to comply with 
established mathematical principles. And so, it is the 
author's contention that successful engineering research fre- 
quently requires a good appreciation of mathematical prin- 
ciples. And these needs will increase, not decrease, as we 
rely more and more on computational tools. 
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